#ifndef __CURLTP23D3D__
#define __CURLTP23D3D__
#include <stdbool.h>
#include "mesh.h"
#include "enumerations.h"
#include "constants.h"
// ***********************************************************************
// Hcurl P2 element, CURL CONFORMING EDGE ELEMENTS OF NÉDÉLEC, 3D
// ***********************************************************************
static INT Curl_T_P2_3D_3D_dof[4] = {0, 2, 2, 0};
static INT Curl_T_P2_3D_3D_Num_Bas = 20;
static INT Curl_T_P2_3D_3D_Value_Dim = 3;
static INT Curl_T_P2_3D_3D_Inter_Dim = 2;
static INT Curl_T_P2_3D_3D_Polydeg = 2;
static bool Curl_T_P2_3D_3D_IsOnlyDependOnRefCoord = 0;
static INT Curl_T_P2_3D_3D_Accuracy = 2;
static MAPTYPE Curl_T_P2_3D_3D_Maptype = Piola;

static void Curl_T_P2_3D_3D_InterFun(DOUBLE RefCoord[], INT dim, DOUBLE *values)
{
   DOUBLE x, y, z;
   x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
   values[0] = 8*x*y - 12*y - 12*z - 6*x + 8*x*z + 16*y*z + 8*y*y + 8*z*z + 4;
   values[1] = 6*x - 8*x*y - 8*x*z - 8*x*x;
   values[2] = 6*x - 8*x*y - 8*x*z - 8*x*x;
   values[3] = 12*x + 14*y + 14*z - 16*x*y - 16*x*z - 16*y*z - 8*y*y - 8*z*z - 6;
   values[4] = 8*x*y - 10*x + 8*x*z + 16*x*x;
   values[5] = 8*x*y - 10*x + 8*x*z + 16*x*x;
   values[6] = 6*y - 8*x*y - 8*y*z - 8*y*y;
   values[7] = 8*x*y - 6*y - 12*z - 12*x + 16*x*z + 8*y*z + 8*x*x + 8*z*z + 4;
   values[8] = 6*y - 8*x*y - 8*y*z - 8*y*y;
   values[9] = 8*x*y - 10*y + 8*y*z + 16*y*y;
   values[10] = 14*x + 12*y + 14*z - 16*x*y - 16*x*z - 16*y*z - 8*x*x - 8*z*z - 6;
   values[11] = 8*x*y - 10*y + 8*y*z + 16*y*y;
   values[12] = 6*z - 8*x*z - 8*y*z - 8*z*z;
   values[13] = 6*z - 8*x*z - 8*y*z - 8*z*z;
   values[14] = 16*x*y - 12*y - 6*z - 12*x + 8*x*z + 8*y*z + 8*x*x + 8*y*y + 4;
   values[15] = 8*x*z - 10*z + 8*y*z + 16*z*z;
   values[16] = 8*x*z - 10*z + 8*y*z + 16*z*z;
   values[17] = 14*x + 14*y + 12*z - 16*x*y - 16*x*z - 16*y*z - 8*x*x - 8*y*y - 6;
   values[18] = 2*y - 8*x*y;
   values[19] = 8*x*x - 4*x;
   values[20] = 0;
   values[21] = 2*y + 8*x*y - 8*y*y;
   values[22] = 2*x + 8*x*y - 8*x*x;
   values[23] = 0;
   values[24] = 2*z - 8*x*z;
   values[25] = 0;
   values[26] = 8*x*x - 4*x;
   values[27] = 2*z + 8*x*z - 8*z*z;
   values[28] = 0;
   values[29] = 2*x + 8*x*z - 8*x*x;
   values[30] = 0;
   values[31] = 2*z - 8*y*z;
   values[32] = 8*y*y - 4*y;
   values[33] = 0;
   values[34] = 2*z + 8*y*z - 8*z*z;
   values[35] = 2*y + 8*y*z - 8*y*y;
   values[36] = -4*y*z;
   values[37] = 8*x*z;
   values[38] = -4*x*y;
   values[39] = -4*y*z;
   values[40] = -4*x*z;
   values[41] = 8*x*y;
   values[42] = 4*y*z;
   values[43] = 8*z - 8*x*z - 4*y*z - 8*z*z;
   values[44] = 4*x*y - 4*y + 8*y*z + 4*y*y;
   values[45] = 4*y*z;
   values[46] = 4*x*z - 4*z + 8*y*z + 4*z*z;
   values[47] = 8*y - 8*x*y - 4*y*z - 8*y*y;
   values[48] = 8*z - 4*x*z - 8*y*z - 8*z*z;
   values[49] = 4*x*z;
   values[50] = 4*x*y - 4*x + 8*x*z + 4*x*x;
   values[51] = 8*x*z - 4*z + 4*y*z + 4*z*z;
   values[52] = 4*x*z;
   values[53] = 8*x - 8*x*y - 4*x*z - 8*x*x;
   values[54] = 8*y - 4*x*y - 8*y*z - 8*y*y;
   values[55] = 8*x*y - 4*x + 4*x*z + 4*x*x;
   values[56] = 4*x*y;
   values[57] = 8*x*y - 4*y + 4*y*z + 4*y*y;
   values[58] = 8*x - 4*x*y - 8*x*z - 8*x*x;
   values[59] = 4*x*y;
}
static void Curl_T_P2_3D_3D_D000(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
   values[0] = 8*x*y - 12*y - 12*z - 6*x + 8*x*z + 16*y*z + 8*y*y + 8*z*z + 4;
   values[1] = 6*x - 8*x*y - 8*x*z - 8*x*x;
   values[2] = 6*x - 8*x*y - 8*x*z - 8*x*x;
   values[3] = 12*x + 14*y + 14*z - 16*x*y - 16*x*z - 16*y*z - 8*y*y - 8*z*z - 6;
   values[4] = 8*x*y - 10*x + 8*x*z + 16*x*x;
   values[5] = 8*x*y - 10*x + 8*x*z + 16*x*x;
   values[6] = 6*y - 8*x*y - 8*y*z - 8*y*y;
   values[7] = 8*x*y - 6*y - 12*z - 12*x + 16*x*z + 8*y*z + 8*x*x + 8*z*z + 4;
   values[8] = 6*y - 8*x*y - 8*y*z - 8*y*y;
   values[9] = 8*x*y - 10*y + 8*y*z + 16*y*y;
   values[10] = 14*x + 12*y + 14*z - 16*x*y - 16*x*z - 16*y*z - 8*x*x - 8*z*z - 6;
   values[11] = 8*x*y - 10*y + 8*y*z + 16*y*y;
   values[12] = 6*z - 8*x*z - 8*y*z - 8*z*z;
   values[13] = 6*z - 8*x*z - 8*y*z - 8*z*z;
   values[14] = 16*x*y - 12*y - 6*z - 12*x + 8*x*z + 8*y*z + 8*x*x + 8*y*y + 4;
   values[15] = 8*x*z - 10*z + 8*y*z + 16*z*z;
   values[16] = 8*x*z - 10*z + 8*y*z + 16*z*z;
   values[17] = 14*x + 14*y + 12*z - 16*x*y - 16*x*z - 16*y*z - 8*x*x - 8*y*y - 6;
   values[18] = 2*y - 8*x*y;
   values[19] = 8*x*x - 4*x;
   values[20] = 0;
   values[21] = 2*y + 8*x*y - 8*y*y;
   values[22] = 2*x + 8*x*y - 8*x*x;
   values[23] = 0;
   values[24] = 2*z - 8*x*z;
   values[25] = 0;
   values[26] = 8*x*x - 4*x;
   values[27] = 2*z + 8*x*z - 8*z*z;
   values[28] = 0;
   values[29] = 2*x + 8*x*z - 8*x*x;
   values[30] = 0;
   values[31] = 2*z - 8*y*z;
   values[32] = 8*y*y - 4*y;
   values[33] = 0;
   values[34] = 2*z + 8*y*z - 8*z*z;
   values[35] = 2*y + 8*y*z - 8*y*y;
   values[36] = -4*y*z;
   values[37] = 8*x*z;
   values[38] = -4*x*y;
   values[39] = -4*y*z;
   values[40] = -4*x*z;
   values[41] = 8*x*y;
   values[42] = 4*y*z;
   values[43] = 8*z - 8*x*z - 4*y*z - 8*z*z;
   values[44] = 4*x*y - 4*y + 8*y*z + 4*y*y;
   values[45] = 4*y*z;
   values[46] = 4*x*z - 4*z + 8*y*z + 4*z*z;
   values[47] = 8*y - 8*x*y - 4*y*z - 8*y*y;
   values[48] = 8*z - 4*x*z - 8*y*z - 8*z*z;
   values[49] = 4*x*z;
   values[50] = 4*x*y - 4*x + 8*x*z + 4*x*x;
   values[51] = 8*x*z - 4*z + 4*y*z + 4*z*z;
   values[52] = 4*x*z;
   values[53] = 8*x - 8*x*y - 4*x*z - 8*x*x;
   values[54] = 8*y - 4*x*y - 8*y*z - 8*y*y;
   values[55] = 8*x*y - 4*x + 4*x*z + 4*x*x;
   values[56] = 4*x*y;
   values[57] = 8*x*y - 4*y + 4*y*z + 4*y*y;
   values[58] = 8*x - 4*x*y - 8*x*z - 8*x*x;
   values[59] = 4*x*y;
}
static void Curl_T_P2_3D_3D_D100(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
   values[0] = 8*y + 8*z - 6;
   values[1] = 6 - 8*y - 8*z - 16*x;
   values[2] = 6 - 8*y - 8*z - 16*x;
   values[3] = 12 - 16*z - 16*y;
   values[4] = 32*x + 8*y + 8*z - 10;
   values[5] = 32*x + 8*y + 8*z - 10;
   values[6] = -8*y;
   values[7] = 16*x + 8*y + 16*z - 12;
   values[8] = -8*y;
   values[9] = 8*y;
   values[10] = 14 - 16*y - 16*z - 16*x;
   values[11] = 8*y;
   values[12] = -8*z;
   values[13] = -8*z;
   values[14] = 16*x + 16*y + 8*z - 12;
   values[15] = 8*z;
   values[16] = 8*z;
   values[17] = 14 - 16*y - 16*z - 16*x;
   values[18] = -8*y;
   values[19] = 16*x - 4;
   values[20] = 0;
   values[21] = 8*y;
   values[22] = 8*y - 16*x + 2;
   values[23] = 0;
   values[24] = -8*z;
   values[25] = 0;
   values[26] = 16*x - 4;
   values[27] = 8*z;
   values[28] = 0;
   values[29] = 8*z - 16*x + 2;
   values[30] = 0;
   values[31] = 0;
   values[32] = 0;
   values[33] = 0;
   values[34] = 0;
   values[35] = 0;
   values[36] = 0;
   values[37] = 8*z;
   values[38] = -4*y;
   values[39] = 0;
   values[40] = -4*z;
   values[41] = 8*y;
   values[42] = 0;
   values[43] = -8*z;
   values[44] = 4*y;
   values[45] = 0;
   values[46] = 4*z;
   values[47] = -8*y;
   values[48] = -4*z;
   values[49] = 4*z;
   values[50] = 8*x + 4*y + 8*z - 4;
   values[51] = 8*z;
   values[52] = 4*z;
   values[53] = 8 - 8*y - 4*z - 16*x;
   values[54] = -4*y;
   values[55] = 8*x + 8*y + 4*z - 4;
   values[56] = 4*y;
   values[57] = 8*y;
   values[58] = 8 - 4*y - 8*z - 16*x;
   values[59] = 4*y;
}
static void Curl_T_P2_3D_3D_D010(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
   values[0] = 8*x + 16*y + 16*z - 12;
   values[1] = -8*x;
   values[2] = -8*x;
   values[3] = 14 - 16*y - 16*z - 16*x;
   values[4] = 8*x;
   values[5] = 8*x;
   values[6] = 6 - 16*y - 8*z - 8*x;
   values[7] = 8*x + 8*z - 6;
   values[8] = 6 - 16*y - 8*z - 8*x;
   values[9] = 8*x + 32*y + 8*z - 10;
   values[10] = 12 - 16*z - 16*x;
   values[11] = 8*x + 32*y + 8*z - 10;
   values[12] = -8*z;
   values[13] = -8*z;
   values[14] = 16*x + 16*y + 8*z - 12;
   values[15] = 8*z;
   values[16] = 8*z;
   values[17] = 14 - 16*y - 16*z - 16*x;
   values[18] = 2 - 8*x;
   values[19] = 0;
   values[20] = 0;
   values[21] = 8*x - 16*y + 2;
   values[22] = 8*x;
   values[23] = 0;
   values[24] = 0;
   values[25] = 0;
   values[26] = 0;
   values[27] = 0;
   values[28] = 0;
   values[29] = 0;
   values[30] = 0;
   values[31] = -8*z;
   values[32] = 16*y - 4;
   values[33] = 0;
   values[34] = 8*z;
   values[35] = 8*z - 16*y + 2;
   values[36] = -4*z;
   values[37] = 0;
   values[38] = -4*x;
   values[39] = -4*z;
   values[40] = 0;
   values[41] = 8*x;
   values[42] = 4*z;
   values[43] = -4*z;
   values[44] = 4*x + 8*y + 8*z - 4;
   values[45] = 4*z;
   values[46] = 8*z;
   values[47] = 8 - 16*y - 4*z - 8*x;
   values[48] = -8*z;
   values[49] = 0;
   values[50] = 4*x;
   values[51] = 4*z;
   values[52] = 0;
   values[53] = -8*x;
   values[54] = 8 - 16*y - 8*z - 4*x;
   values[55] = 8*x;
   values[56] = 4*x;
   values[57] = 8*x + 8*y + 4*z - 4;
   values[58] = -4*x;
   values[59] = 4*x;
}
static void Curl_T_P2_3D_3D_D001(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
   values[0] = 8*x + 16*y + 16*z - 12;
   values[1] = -8*x;
   values[2] = -8*x;
   values[3] = 14 - 16*y - 16*z - 16*x;
   values[4] = 8*x;
   values[5] = 8*x;
   values[6] = -8*y;
   values[7] = 16*x + 8*y + 16*z - 12;
   values[8] = -8*y;
   values[9] = 8*y;
   values[10] = 14 - 16*y - 16*z - 16*x;
   values[11] = 8*y;
   values[12] = 6 - 8*y - 16*z - 8*x;
   values[13] = 6 - 8*y - 16*z - 8*x;
   values[14] = 8*x + 8*y - 6;
   values[15] = 8*x + 8*y + 32*z - 10;
   values[16] = 8*x + 8*y + 32*z - 10;
   values[17] = 12 - 16*y - 16*x;
   values[18] = 0;
   values[19] = 0;
   values[20] = 0;
   values[21] = 0;
   values[22] = 0;
   values[23] = 0;
   values[24] = 2 - 8*x;
   values[25] = 0;
   values[26] = 0;
   values[27] = 8*x - 16*z + 2;
   values[28] = 0;
   values[29] = 8*x;
   values[30] = 0;
   values[31] = 2 - 8*y;
   values[32] = 0;
   values[33] = 0;
   values[34] = 8*y - 16*z + 2;
   values[35] = 8*y;
   values[36] = -4*y;
   values[37] = 8*x;
   values[38] = 0;
   values[39] = -4*y;
   values[40] = -4*x;
   values[41] = 0;
   values[42] = 4*y;
   values[43] = 8 - 4*y - 16*z - 8*x;
   values[44] = 8*y;
   values[45] = 4*y;
   values[46] = 4*x + 8*y + 8*z - 4;
   values[47] = -4*y;
   values[48] = 8 - 8*y - 16*z - 4*x;
   values[49] = 4*x;
   values[50] = 8*x;
   values[51] = 8*x + 4*y + 8*z - 4;
   values[52] = 4*x;
   values[53] = -4*x;
   values[54] = -8*y;
   values[55] = 4*x;
   values[56] = 0;
   values[57] = 4*y;
   values[58] = -8*x;
   values[59] = 0;
}
static void Curl_T_P2_3D_3D_D200(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
   values[0] = 0;
   values[1] = -16;
   values[2] = -16;
   values[3] = 0;
   values[4] = 32;
   values[5] = 32;
   values[6] = 0;
   values[7] = 16;
   values[8] = 0;
   values[9] = 0;
   values[10] = -16;
   values[11] = 0;
   values[12] = 0;
   values[13] = 0;
   values[14] = 16;
   values[15] = 0;
   values[16] = 0;
   values[17] = -16;
   values[18] = 0;
   values[19] = 16;
   values[20] = 0;
   values[21] = 0;
   values[22] = -16;
   values[23] = 0;
   values[24] = 0;
   values[25] = 0;
   values[26] = 16;
   values[27] = 0;
   values[28] = 0;
   values[29] = -16;
   values[30] = 0;
   values[31] = 0;
   values[32] = 0;
   values[33] = 0;
   values[34] = 0;
   values[35] = 0;
   values[36] = 0;
   values[37] = 0;
   values[38] = 0;
   values[39] = 0;
   values[40] = 0;
   values[41] = 0;
   values[42] = 0;
   values[43] = 0;
   values[44] = 0;
   values[45] = 0;
   values[46] = 0;
   values[47] = 0;
   values[48] = 0;
   values[49] = 0;
   values[50] = 8;
   values[51] = 0;
   values[52] = 0;
   values[53] = -16;
   values[54] = 0;
   values[55] = 8;
   values[56] = 0;
   values[57] = 0;
   values[58] = -16;
   values[59] = 0;
}
static void Curl_T_P2_3D_3D_D020(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
   values[0] = 16;
   values[1] = 0;
   values[2] = 0;
   values[3] = -16;
   values[4] = 0;
   values[5] = 0;
   values[6] = -16;
   values[7] = 0;
   values[8] = -16;
   values[9] = 32;
   values[10] = 0;
   values[11] = 32;
   values[12] = 0;
   values[13] = 0;
   values[14] = 16;
   values[15] = 0;
   values[16] = 0;
   values[17] = -16;
   values[18] = 0;
   values[19] = 0;
   values[20] = 0;
   values[21] = -16;
   values[22] = 0;
   values[23] = 0;
   values[24] = 0;
   values[25] = 0;
   values[26] = 0;
   values[27] = 0;
   values[28] = 0;
   values[29] = 0;
   values[30] = 0;
   values[31] = 0;
   values[32] = 16;
   values[33] = 0;
   values[34] = 0;
   values[35] = -16;
   values[36] = 0;
   values[37] = 0;
   values[38] = 0;
   values[39] = 0;
   values[40] = 0;
   values[41] = 0;
   values[42] = 0;
   values[43] = 0;
   values[44] = 8;
   values[45] = 0;
   values[46] = 0;
   values[47] = -16;
   values[48] = 0;
   values[49] = 0;
   values[50] = 0;
   values[51] = 0;
   values[52] = 0;
   values[53] = 0;
   values[54] = -16;
   values[55] = 0;
   values[56] = 0;
   values[57] = 8;
   values[58] = 0;
   values[59] = 0;
}
static void Curl_T_P2_3D_3D_D002(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
   values[0] = 16;
   values[1] = 0;
   values[2] = 0;
   values[3] = -16;
   values[4] = 0;
   values[5] = 0;
   values[6] = 0;
   values[7] = 16;
   values[8] = 0;
   values[9] = 0;
   values[10] = -16;
   values[11] = 0;
   values[12] = -16;
   values[13] = -16;
   values[14] = 0;
   values[15] = 32;
   values[16] = 32;
   values[17] = 0;
   values[18] = 0;
   values[19] = 0;
   values[20] = 0;
   values[21] = 0;
   values[22] = 0;
   values[23] = 0;
   values[24] = 0;
   values[25] = 0;
   values[26] = 0;
   values[27] = -16;
   values[28] = 0;
   values[29] = 0;
   values[30] = 0;
   values[31] = 0;
   values[32] = 0;
   values[33] = 0;
   values[34] = -16;
   values[35] = 0;
   values[36] = 0;
   values[37] = 0;
   values[38] = 0;
   values[39] = 0;
   values[40] = 0;
   values[41] = 0;
   values[42] = 0;
   values[43] = -16;
   values[44] = 0;
   values[45] = 0;
   values[46] = 8;
   values[47] = 0;
   values[48] = -16;
   values[49] = 0;
   values[50] = 0;
   values[51] = 8;
   values[52] = 0;
   values[53] = 0;
   values[54] = 0;
   values[55] = 0;
   values[56] = 0;
   values[57] = 0;
   values[58] = 0;
   values[59] = 0;
}
static void Curl_T_P2_3D_3D_D110(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
   values[0] = 8;
   values[1] = -8;
   values[2] = -8;
   values[3] = -16;
   values[4] = 8;
   values[5] = 8;
   values[6] = -8;
   values[7] = 8;
   values[8] = -8;
   values[9] = 8;
   values[10] = -16;
   values[11] = 8;
   values[12] = 0;
   values[13] = 0;
   values[14] = 16;
   values[15] = 0;
   values[16] = 0;
   values[17] = -16;
   values[18] = -8;
   values[19] = 0;
   values[20] = 0;
   values[21] = 8;
   values[22] = 8;
   values[23] = 0;
   values[24] = 0;
   values[25] = 0;
   values[26] = 0;
   values[27] = 0;
   values[28] = 0;
   values[29] = 0;
   values[30] = 0;
   values[31] = 0;
   values[32] = 0;
   values[33] = 0;
   values[34] = 0;
   values[35] = 0;
   values[36] = 0;
   values[37] = 0;
   values[38] = -4;
   values[39] = 0;
   values[40] = 0;
   values[41] = 8;
   values[42] = 0;
   values[43] = 0;
   values[44] = 4;
   values[45] = 0;
   values[46] = 0;
   values[47] = -8;
   values[48] = 0;
   values[49] = 0;
   values[50] = 4;
   values[51] = 0;
   values[52] = 0;
   values[53] = -8;
   values[54] = -4;
   values[55] = 8;
   values[56] = 4;
   values[57] = 8;
   values[58] = -4;
   values[59] = 4;
}
static void Curl_T_P2_3D_3D_D101(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
   values[0] = 8;
   values[1] = -8;
   values[2] = -8;
   values[3] = -16;
   values[4] = 8;
   values[5] = 8;
   values[6] = 0;
   values[7] = 16;
   values[8] = 0;
   values[9] = 0;
   values[10] = -16;
   values[11] = 0;
   values[12] = -8;
   values[13] = -8;
   values[14] = 8;
   values[15] = 8;
   values[16] = 8;
   values[17] = -16;
   values[18] = 0;
   values[19] = 0;
   values[20] = 0;
   values[21] = 0;
   values[22] = 0;
   values[23] = 0;
   values[24] = -8;
   values[25] = 0;
   values[26] = 0;
   values[27] = 8;
   values[28] = 0;
   values[29] = 8;
   values[30] = 0;
   values[31] = 0;
   values[32] = 0;
   values[33] = 0;
   values[34] = 0;
   values[35] = 0;
   values[36] = 0;
   values[37] = 8;
   values[38] = 0;
   values[39] = 0;
   values[40] = -4;
   values[41] = 0;
   values[42] = 0;
   values[43] = -8;
   values[44] = 0;
   values[45] = 0;
   values[46] = 4;
   values[47] = 0;
   values[48] = -4;
   values[49] = 4;
   values[50] = 8;
   values[51] = 8;
   values[52] = 4;
   values[53] = -4;
   values[54] = 0;
   values[55] = 4;
   values[56] = 0;
   values[57] = 0;
   values[58] = -8;
   values[59] = 0;
}
static void Curl_T_P2_3D_3D_D011(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y, z;
   x=RefCoord[0]; y=RefCoord[1]; z=RefCoord[2];
   values[0] = 16;
   values[1] = 0;
   values[2] = 0;
   values[3] = -16;
   values[4] = 0;
   values[5] = 0;
   values[6] = -8;
   values[7] = 8;
   values[8] = -8;
   values[9] = 8;
   values[10] = -16;
   values[11] = 8;
   values[12] = -8;
   values[13] = -8;
   values[14] = 8;
   values[15] = 8;
   values[16] = 8;
   values[17] = -16;
   values[18] = 0;
   values[19] = 0;
   values[20] = 0;
   values[21] = 0;
   values[22] = 0;
   values[23] = 0;
   values[24] = 0;
   values[25] = 0;
   values[26] = 0;
   values[27] = 0;
   values[28] = 0;
   values[29] = 0;
   values[30] = 0;
   values[31] = -8;
   values[32] = 0;
   values[33] = 0;
   values[34] = 8;
   values[35] = 8;
   values[36] = -4;
   values[37] = 0;
   values[38] = 0;
   values[39] = -4;
   values[40] = 0;
   values[41] = 0;
   values[42] = 4;
   values[43] = -4;
   values[44] = 8;
   values[45] = 4;
   values[46] = 8;
   values[47] = -4;
   values[48] = -8;
   values[49] = 0;
   values[50] = 0;
   values[51] = 4;
   values[52] = 0;
   values[53] = 0;
   values[54] = -8;
   values[55] = 0;
   values[56] = 0;
   values[57] = 4;
   values[58] = 0;
   values[59] = 0;
}

static void Curl_T_P2_3D_3D_Nodal(ELEMENT *elem, FUNCTIONVEC *fun, INT dim, DOUBLE *values)
{
   // 直接给出积分格式
   INT NumPoints = 4;
   DOUBLE QuadX14[4] = {0.069431844202973713731097404888715,
                        0.33000947820757187134432797392947,
                        0.66999052179242812865567202607053,
                        0.93056815579702628626890259511129};
   DOUBLE QuadW14[4] = {0.17392742256872692485636378,
                        0.3260725774312730473880606,
                        0.3260725774312730473880606,
                        0.17392742256872692485636378};
   INT vert0[6] = {0, 0, 0, 1, 1, 2};
   INT vert1[6] = {1, 2, 3, 2, 3, 3};
   INT i, j, k;
   INT num_fun = dim / 3; // 输入函数的个数(curl元是3D的)
   DOUBLE coord[3];
   DOUBLE vec_line[3]; // 带方向的线(按照单元内部的方向)
   DOUBLE *values_temp = malloc(dim * sizeof(DOUBLE));
#if NEDELECSCALE == 0
   for (i = 0; i < 6; i++)
   {
      for (k = 0; k < num_fun; k++)
      {
         values[k + 2 * i * num_fun] = 0.0;
         values[k + (2 * i + 1) * num_fun] = 0.0;
      }
      vec_line[0] = elem->Vert_X[vert1[i]] - elem->Vert_X[vert0[i]];
      vec_line[1] = elem->Vert_Y[vert1[i]] - elem->Vert_Y[vert0[i]];
      vec_line[2] = elem->Vert_Z[vert1[i]] - elem->Vert_Z[vert0[i]];
      for (j = 0; j < NumPoints; j++)
      {
         // coord
         coord[0] = (1 - QuadX14[j]) * elem->Vert_X[vert0[i]] + QuadX14[j] * elem->Vert_X[vert1[i]];
         coord[1] = (1 - QuadX14[j]) * elem->Vert_Y[vert0[i]] + QuadX14[j] * elem->Vert_Y[vert1[i]];
         coord[2] = (1 - QuadX14[j]) * elem->Vert_Z[vert0[i]] + QuadX14[j] * elem->Vert_Z[vert1[i]];
         // u(t_j)
         fun(coord, dim, values_temp);
         // 线上有两个自由度
         // \int_{e} \vec{u}\cdot \vec{\tau}de = \sum_j \vec{u}(t_j)cdot \vec{\tau}\omega_j|e| = \sum_j \vec{u}(t_j)cdot vec_line \omega_j
         // \int_{e} \vec{u}\cdot \vec{\tau}de = \sum_j \vec{u}(t_j)cdot \vec{\tau} t \omega_j|e| = \sum_j \vec{u}(t_j)cdot vec_line \omega_j \X_j
         for (k = 0; k < num_fun; k++)
         {
            values[k + 2 * i * num_fun] += (values_temp[3 * k] * vec_line[0] + values_temp[3 * k + 1] * vec_line[1] + values_temp[3 * k + 2] * vec_line[2]) * QuadW14[j];
            values[k + (2 * i + 1) * num_fun] += (values_temp[3 * k] * vec_line[0] + values_temp[3 * k + 1] * vec_line[1] + values_temp[3 * k + 2] * vec_line[2]) * QuadW14[j] * QuadX14[j];
         }
      }
   }
#else
   DOUBLE line_length[6];
   for (i = 0; i < 6; i++)
   {
      for (k = 0; k < num_fun; k++)
      {
         values[k + 2 * i * num_fun] = 0.0;
         values[k + (2 * i + 1) * num_fun] = 0.0;
      }
      vec_line[0] = elem->Vert_X[vert1[i]] - elem->Vert_X[vert0[i]];
      vec_line[1] = elem->Vert_Y[vert1[i]] - elem->Vert_Y[vert0[i]];
      vec_line[2] = elem->Vert_Z[vert1[i]] - elem->Vert_Z[vert0[i]];
      line_length[i] = sqrt(vec_line[0]*vec_line[0]+vec_line[1]*vec_line[1]+vec_line[2]*vec_line[2]);
      for (j = 0; j < NumPoints; j++)
      {
         // coord
         coord[0] = (1 - QuadX14[j]) * elem->Vert_X[vert0[i]] + QuadX14[j] * elem->Vert_X[vert1[i]];
         coord[1] = (1 - QuadX14[j]) * elem->Vert_Y[vert0[i]] + QuadX14[j] * elem->Vert_Y[vert1[i]];
         coord[2] = (1 - QuadX14[j]) * elem->Vert_Z[vert0[i]] + QuadX14[j] * elem->Vert_Z[vert1[i]];
         // u(t_j)
         fun(coord, dim, values_temp);
         // 线上有两个自由度
         // \int_{e} \vec{u}\cdot \vec{\tau}de = \sum_j \vec{u}(t_j)cdot \vec{\tau}\omega_j|e| = \sum_j \vec{u}(t_j)cdot vec_line \omega_j
         // \int_{e} \vec{u}\cdot \vec{\tau}de = \sum_j \vec{u}(t_j)cdot \vec{\tau} t \omega_j|e| = \sum_j \vec{u}(t_j)cdot vec_line \omega_j \X_j
         for (k = 0; k < num_fun; k++)
         {
            values[k + 2 * i * num_fun] += (values_temp[3 * k] * vec_line[0] + values_temp[3 * k + 1] * vec_line[1] + values_temp[3 * k + 2] * vec_line[2]) * QuadW14[j] / line_length[i];
            values[k + (2 * i + 1) * num_fun] += (values_temp[3 * k] * vec_line[0] + values_temp[3 * k + 1] * vec_line[1] + values_temp[3 * k + 2] * vec_line[2]) * QuadW14[j] * QuadX14[j] / line_length[i];
         }
      }
   }
#endif
   NumPoints = 13;
   DOUBLE QuadX213[13] = {0.333333333333000,
                          0.260345966079000,
                          0.260345966079000,
                          0.479308067842000,
                          0.065130102902000,
                          0.065130102902000,
                          0.869739794196000,
                          0.312865496005000,
                          0.638444188570000,
                          0.048690315425000,
                          0.638444188570000,
                          0.312865496005000,
                          0.048690315425000};
   DOUBLE QuadY213[13] = {0.333333333333000,
                          0.260345966079000,
                          0.479308067842000,
                          0.260345966079000,
                          0.065130102902000,
                          0.869739794196000,
                          0.065130102902000,
                          0.638444188570000,
                          0.048690315425000,
                          0.312865496005000,
                          0.312865496005000,
                          0.048690315425000,
                          0.638444188570000};
   DOUBLE QuadW213[13] = {-0.149570044468000,
                          0.175615257433000,
                          0.175615257433000,
                          0.175615257433000,
                          0.053347235609000,
                          0.053347235609000,
                          0.053347235609000,
                          0.077113760890000,
                          0.077113760890000,
                          0.077113760890000,
                          0.077113760890000,
                          0.077113760890000,
                          0.077113760890000};
   INT vert04face[4] = {1, 0, 0, 0};
   INT vert14face[4] = {2, 2, 1, 1};
   INT vert24face[4] = {3, 3, 3, 2};
   DOUBLE tau0[3]; // 面的第一条边的方向
   DOUBLE tau1[3]; // 面的第二条边的方向
   DOUBLE length[3];
   DOUBLE area;
#if NEDELECSCALE == 0
   for (i = 0; i < 4; i++)
   {
      for (k = 0; k < num_fun; k++)
      {
         values[k + 2 * (i + 6) * num_fun] = 0.0;
         values[k + (2 * (i + 6) + 1) * num_fun] = 0.0;
      }
      // 计算面的第一条边的方向
      tau0[0] = elem->Vert_X[vert14face[i]] - elem->Vert_X[vert04face[i]];
      tau0[1] = elem->Vert_Y[vert14face[i]] - elem->Vert_Y[vert04face[i]];
      tau0[2] = elem->Vert_Z[vert14face[i]] - elem->Vert_Z[vert04face[i]];
      // 计算面的第二条边的方向
      tau1[0] = elem->Vert_X[vert24face[i]] - elem->Vert_X[vert04face[i]];
      tau1[1] = elem->Vert_Y[vert24face[i]] - elem->Vert_Y[vert04face[i]];
      tau1[2] = elem->Vert_Z[vert24face[i]] - elem->Vert_Z[vert04face[i]];
      for (j = 0; j < NumPoints; j++)
      {
         // coord
         coord[0] = (1.0 - QuadX213[j] - QuadY213[j]) * elem->Vert_X[vert04face[i]] + QuadX213[j] * elem->Vert_X[vert14face[i]] + QuadY213[j] * elem->Vert_X[vert24face[i]];
         coord[1] = (1.0 - QuadX213[j] - QuadY213[j]) * elem->Vert_Y[vert04face[i]] + QuadX213[j] * elem->Vert_Y[vert14face[i]] + QuadY213[j] * elem->Vert_Y[vert24face[i]];
         coord[2] = (1.0 - QuadX213[j] - QuadY213[j]) * elem->Vert_Z[vert04face[i]] + QuadX213[j] * elem->Vert_Z[vert14face[i]] + QuadY213[j] * elem->Vert_Z[vert24face[i]];
         // u(t_j)
         fun(coord, dim, values_temp);
         // 面上有两个自由度
         // \int_{f} \vec{u}\cdot \vec{\tau_0}df = \sum_j \vec{u}(t_j)cdot \vec{\tau_0}\omega_j|f| = \sum_j \vec{u}(t_j)cdot \vec{\tau_0} \omega_j|f|
         // \int_{f} \vec{u}\cdot \vec{\tau_1}df = \sum_j \vec{u}(t_j)cdot \vec{\tau_1}\omega_j|f| = \sum_j \vec{u}(t_j)cdot \vec{\tau_1} \omega_j|f|
         for (k = 0; k < num_fun; k++)
         {
            values[k + 2 * (i + 6) * num_fun] += (values_temp[3 * k] * tau0[0] + values_temp[3 * k + 1] * tau0[1] + values_temp[3 * k + 2] * tau0[2]) * QuadW213[j];
            values[k + (2 * (i + 6) + 1) * num_fun] += (values_temp[3 * k] * tau1[0] + values_temp[3 * k + 1] * tau1[1] + values_temp[3 * k + 2] * tau1[2]) * QuadW213[j];
         }
      }
   }
#else
   INT line4face1[4] = {3,1,0,0};
   INT line4face2[4] = {4,2,2,1};
   for (i = 0; i < 4; i++)
   {
      for (k = 0; k < num_fun; k++)
      {
         values[k + 2 * (i + 6) * num_fun] = 0.0;
         values[k + (2 * (i + 6) + 1) * num_fun] = 0.0;
      }
      // 计算面的第一条边的方向
      tau0[0] = elem->Vert_X[vert14face[i]] - elem->Vert_X[vert04face[i]];
      tau0[1] = elem->Vert_Y[vert14face[i]] - elem->Vert_Y[vert04face[i]];
      tau0[2] = elem->Vert_Z[vert14face[i]] - elem->Vert_Z[vert04face[i]];
      // 计算面的第二条边的方向
      tau1[0] = elem->Vert_X[vert24face[i]] - elem->Vert_X[vert04face[i]];
      tau1[1] = elem->Vert_Y[vert24face[i]] - elem->Vert_Y[vert04face[i]];
      tau1[2] = elem->Vert_Z[vert24face[i]] - elem->Vert_Z[vert04face[i]];
      for (j = 0; j < NumPoints; j++)
      {
         // coord
         coord[0] = (1.0 - QuadX213[j] - QuadY213[j]) * elem->Vert_X[vert04face[i]] + QuadX213[j] * elem->Vert_X[vert14face[i]] + QuadY213[j] * elem->Vert_X[vert24face[i]];
         coord[1] = (1.0 - QuadX213[j] - QuadY213[j]) * elem->Vert_Y[vert04face[i]] + QuadX213[j] * elem->Vert_Y[vert14face[i]] + QuadY213[j] * elem->Vert_Y[vert24face[i]];
         coord[2] = (1.0 - QuadX213[j] - QuadY213[j]) * elem->Vert_Z[vert04face[i]] + QuadX213[j] * elem->Vert_Z[vert14face[i]] + QuadY213[j] * elem->Vert_Z[vert24face[i]];
         // u(t_j)
         fun(coord, dim, values_temp);
         // 面上有两个自由度
         // \int_{f} \vec{u}\cdot \vec{\tau_0}df = \sum_j \vec{u}(t_j)cdot \vec{\tau_0}\omega_j|f| = \sum_j \vec{u}(t_j)cdot \vec{\tau_0} \omega_j|f|
         // \int_{f} \vec{u}\cdot \vec{\tau_1}df = \sum_j \vec{u}(t_j)cdot \vec{\tau_1}\omega_j|f| = \sum_j \vec{u}(t_j)cdot \vec{\tau_1} \omega_j|f|
         for (k = 0; k < num_fun; k++)
         {
            values[k + 2 * (i + 6) * num_fun] += (values_temp[3 * k] * tau0[0] + values_temp[3 * k + 1] * tau0[1] + values_temp[3 * k + 2] * tau0[2]) * QuadW213[j] / line_length[line4face1[i]];
            values[k + (2 * (i + 6) + 1) * num_fun] += (values_temp[3 * k] * tau1[0] + values_temp[3 * k + 1] * tau1[1] + values_temp[3 * k + 2] * tau1[2]) * QuadW213[j] / line_length[line4face2[i]];
         }
      }
   }
#endif
   free(values_temp);
}
#endif
